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1.      INTRODUCTION 

The  method  for  solving  systems    of  nonlinear  algebraic  equations 
that   is  the  subject   of  this   investigation  was    developed  by  van  Me  lie    [l]. 
His   subroutine  DIFMF3  was    designed  to  solve  the  steady-state  problem  of 
ordinary   differential  equations  ,   and  for  this   purpose   each   component  of 
the   vector  to  be   solved  for   could  be   either   a  derivative   or   an   ordinary 
variable.      The    development   of  this   method  grew  out   of  ideas  presented  by 
Gear    [2,3]    concerning  the   solution   of  systems    of  mixed  algebraic   and 
differential  equations. 

In   order  to  evaluate   this   method  simply   as    a  means   of  solving 
nonlinear  equations  ,    another  subroutine    called  INTECH  was  prepared  by 
the    author.      This   subroutine   is    essentially   the   same   as    DIFMF3  but   is 
simplified  so   as   to  handle   only  a  single  vector  of  variables   instead  of 
giving  the   option   of  choosing  each   component   from  either  a  vector   of 
variables   or  a  vector  of   derivatives.      Provisions    for  applying  constraints 
to  the   variables  were   also   omitted,    along  with   certain   variable  names   that 
were   used  in   DIFMF3  only  to  pass    data  to  other  subroutines   in   the 
differential -equation-solving  package   of  which   it  was    a  part.      Finally, 
certain    changes  were   required  to   convert   the  subroutine   from  FORTRAN   IV 
to  CALL  OS   FORTRAN,  which  the   author  chose  to   use  because   of  the 
convenience   of  the   timesharing  system  available   for  this    language. 

Most   of  the  examples   used  in  the  investigation  were  taken   from 
articles  by  Boggs    [U,5]   or  by  Powell    [6,7]   so  that   comparisons    could  be 
made  with  their  published  results.      Their  results    appeared  to  the   author 
to  be   good  ones   to    compare  with  because  they  both   report  on   recent 
developments,   their  methods    are   different,    and  both  show  some   improvement 


. _„__!  ■'_.-"* 


2.      SUBROUTINE   INTECH 

This    section   presents    a   description   of  the  "basic  method  employed, 
the   strategy   used  in    an    attempt   to  economize    on  the    amount   of   computation 
required  to   obtain    a  solution,    and  a  discussion   of  some   of  the   details    to 
be   considered  in  using  the  subroutine.      For  a  listing  of  the   subroutine, 
along  with  the    main  program  and  subroutines    used  for  one    of  the   examples, 
refer  to   the  Appendix. 

2. 1     Description    of  the   Method 

The   problem  that  we  wish   to  solve  may  be  written   in  the   form 

f(y)   =  0  (2.1) 

r    1       2  niT 

where  y  =    Ly    ,  y    >•  •  ■   y    J 

12  _n  T 

and  f(y)    =    [f   (y),    f   (y),...    i    (y)]      are  n-component   column  vectors,    and 

0   is    an  n-component   column   vector  of   zeros.      Also,  we  will  use  J(y)    to 

denote   the   Jacobian   matrix 

J(y)  -&&k  . 

y 

Following  Boggs '    approach   in  Section  2   of  his   paper   [U],  we   consider  the 
operator  homotopy  G(s,y)    =  0  which   imbeds    the   real  parameter  s   into  the 
original  equation   so  that   G(0,y)    =  0  has  the  solution  y    ,   and 
G(s,y)     ►   f(y)    as   s  ■*-  °°.      Specifically,  we   use 

G(s,y)    =    f(y)    -  e"Sf(y0)    =  0 
to  achieve   these   results.      Then,   since   G(s,y)   =  0,  we   may  obtain  y' 


dy_ 
ds 

rale ,    in  the   form 


(denoting  r~-  )   by   differentiating  with   respect   to   s   and  using  the    chain 


y'    =  -J   1(y)f(y),  y(0)   =  yQ.  (2.2) 


Euler's   method  of  numerical  integration   applied  to  this    differential 
equation   gives  the  iteration 

ym  =  ym-i +  hy;-i  (2-3) 

We  note   that  this   is   equivalent,  if  we   make  h  =   1,   to  Newton's   method 
applied  to  the   original  equation   (2.1).      Formula  (2.3)    is   used  as   the 
predictor  step   in  van  Melle's   method.      For  the   corrector,  he  uses 

ym  =  ym-l+h(aym+    (l"a)    K-l]  {2-k) 

with   0   <_  a  £  1.     When  a  =  0,    (2.U)    is  the  same  as    (2.3).      If  a  =   0 

and  h  =   1,   the   iteration   of  the  predictor   and  corrector  steps   is   identical 

to  Newton's   method  applied  to    (2.1),   and  if  it    converges  to   a  solution, 

will  do   so   relatively  quickly.      On  the   other  hand,   if  a  =   1  and  h   <<  1, 

the  method  is  no  longer  Newton's   method  and  will  generally   require   many 

steps   to   converge,  but  the    convergence   is  more   likely  to   occur. 

Note  that,    as   it  stands,    (2.1+)    cannot  be   solved  directly   for 

y     because  y',  which  is  a  function  of  y     by   (2.2),   appears   on  the  right- 
m  m  m 

hand  side.      Rather  than  solving  (2.k)    iteratively  each  time,   thus    getting 

an  implicit   method,   we   obtain   an   approximate   solution   in   an  explicit 

form.      To   do  this,    denote  the   right-hand  side   of   (2.3),  which  is  the 

result  of  the  predictor  operation,  by  y        ,   then   denote  the   right-hand 

side   of   (2.1+)   by  y     =  y  +  Ay    .      Using  a  first    order  Taylor  expansion 

mm  m 

on   the   right-hand  side   of   (2.2)    to   represent  y' ,   the   following  explicit 
corrector  may  be  obtained  so  that  the  method  becomes 


A  ■  est  thJ_1  <0>> +  ^-J  <2-6» 

(O)  ,  V 

ym=ym      -  «a  (2.7) 

hy;  =  hy;.i-A  (2'8) 

where    (2.5)    is  the  predictor   and   (2.6)   through    (2.8)    make  up  the 

corrector.      A  single   iteration   of  the   corrector  is   used  after  each 

application   of  the   predictor. 

If  the   inverse  Jacobian  were   to  "be   evaluated  at  every  step, 

we  would  use   J      (y        )    in    (2.6).      However,  part  of  the  strategy  in 
m 

applying  this   method  is   an   attempt   to  reduce   the   amount  of  computation 
required  by  evaluating  the  Jacobian   and  its   inverse   infrequently,    as 
discussed  in  Section  2.2.      Therefore,   the  symbol  J        is   used  in    (2.6) 
to  denote   the  inverse  Jacobian  that   is    currently  being  used. 


2.2      Strategy 

One   of  the   major  problems   that   had  to  be   solved  empirically 
was   the  method   of  varying  a   and  h   in  order  to   find  a  solution   as 
efficiently   as   possible.      The  method  used  in   INTECH   is  basically  the 
same   as  that   used  in  van   Me  lie's   subroutine   DIFMF3    [l],  with   any 
differences   being  noted  in  the   description   that   follows.      In   order  to 
relate  the   discussion  more   closely  to  the   listing  of  the   subroutine 
INTECH   in  the  Appendix,  the  variable  names   ALPHA  and  H  will  be   used  to 
denote,    respectively,   the  symbols   a  and  h  used  in  Section  2.1. 


Initially,   the  procedure  begins  with  ALPHA  =  0   and  H  =   1, 
so  that   the   method  is    the  same   as   Newton's   method.      As   long  as  the  norm 

l|f(y)M  =   *    ly1! 

i=l 
decreases   significantly,   ALPHA  and  H  are  not   changed.      The  norm  is 
calculated  following  each  predictor  step.      If  the  norm  fails   to   decrease 
by   at    least   five  percent   from  its   previous   value,   the   replacements 
ALPHA  =  1   and  H  =    .01   are  made   and  the  previous  vector  y  is    restored 
before   repeating  the   step.      After  this,   on  each  succeeding  iteration   as 
long  as   the  norm  has    decreased  by   at  least  two  percent   from  its  previous 
value,  ALPHA  is    decreased  by   multiplying  it  by  the   factor  0.8,    and  H  is 
increased  by  multiplying  it  by  the   factor  R  =   1.7   -   .  85H  +    .  15/H  which 
is    designed  so  that  the   limitation  H  <_  1  will  be   observed.      This   is   one 
point  at  which  INTECH   differs   from  van  Me  lie's    subroutine  DIFMF3.      In 
DIFMF3  the   corresponding  formula  for  R  is:      R  =    .6  +    .k/E.      In  preliminary 
tests    on  INTECH  involving  just  Examples    1  through  9»   it  was    found  that 
when   using  the   latter  formula  for  R,  no    convergence  was    obtained  for 
Examples    5   and  9-      The  indicated  formula  for  R  was    found  to  be   the  best 
of  several   that  were  tried,    giving  convergence   for  all  nine  examples    and 
requiring  the  same    or   a  smaller  number  of  function  evaluations  than  the 
original  formula  for   all  except  Example    39   in  which   case   one   additional 
evaluation  was   required. 

If  on  any  iteration  after  the  Initially  used  Newton's  method 
has    failed,   the  norm  exceeds  the  previous   norm  by   a  factor  of  100   or  more, 
the   replacements  ALPHA  =   1  and  H  =  max(H/2,   0.2)    are  made,   the  previous 
vector  y  is   restored,   and  the  predictor  step  is    repeated.      However,   if  the 
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ratio  of  the  norm  to  the  previous  norm  remains  between   O.98  and  100, 
the    replacements   ALPHA  =   1   and  R  =  min(l.3,    . 6/H)    are   made    and  the 
corrector  step  is  performed. 

Another  feature   of  the  program  related  to   strategy  is   the 
separation    of  the  vector  y   of  variables    into   two   vectors,    given   the 
names    Y   and  YL  in   the   subroutine   INTECH.      The   values    of   all  variables 
that   appear  nonlinear ly  in   one    or  more   of  the   functions   are   stored  in 
Y;  while   those   that   appear  only  linearly   are  stored  in  YL.      This 
separation  makes   it    convenient  to  take   advantage  of  the   fact,  which 
is   discussed  by  Gear    [2,3],   that   only  the    correction  process    and  not 
the   prediction  process   needs   to  be   applied  to  the   linear  variables. 
For  systems    involving  a  large  number  of  linear  variables  ,  this    feature 
may   reduce  the   amount    of   computation   considerably. 

Another  problem  that  had  to  be   considered  was   the   frequency  of 
evaluating  and  inverting  the  Jacobian.      It   is   not    desirable   to   do   this 
on  every  step,   particularly  where  n   is    large,   because   of  the  excessive 
amount  of  computation  required.      The  strategy  that  has  been   developed 
calls    for  recomputing  the  Jacobian  every   5NY   steps  where  NY   is  the 
number  of  variables   in   the  vector  Y   discussed  in   the  preceding  paragraph. 
However,    if  the  Jacobian  has  not  been  evaluated  recently   and,   in 
addition,   either  the  norm  increases   by  a  factor  of  100   or  more   from  its 
previous   value   or  else   it  becomes   less   than   1.0    (indicating  that  a 
solution   is  being  approached)  ,  then   the  Jacobian   is   computed  and  inverted 
even  though  this   action  is  not   called  for  by   the   first    rule. 
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2  .3     The  Calling  Program  and  Other  Subroutines 

The   calling  program  must    define   the   variables   required  by  the 
parameter  list,  which  is   given   in  the   listing  of  INTECH  in  the  Appendix. 
Appropriate    input   values   must  be   supplied  for  NSEND,   the  maximum  number 
of  iterations  to  be  performed;   DEL,   the  value  used  to   compare  with  the 
norm 

n 

*    |yxl 

i=l 
to   determine  whether  or  not   satisfactory   convergence  has  been   obtained; 
NY,  the  number  of  nonlinear  variables    (the   dimension   of  Y);  NL ,   the 
number  of  linear  variables    (the    dimension   of  YL);   N,   the  total  number 
of  variables,   NY  +  NL ;  NY  starting  values   in  the   array  Y,    and  NL 
starting  values   in  the   array  YL.      Results   are  printed  by  statements 
within  INTECH   after  each  prediction  process   of  the  iteration,   so  no 
provisions    for   output   are  normally  needed  in  the   calling  program.      The 
formats    used  in  INTECH  were    designed  for  use  with   a  terminal  taking 
paper  8  1/2"  wide,    and  should  be   changed  appropriately  if  a  line 
printer  is   to  be   used. 

In   addition  to  the   calling  program,   the   user  must   supply 
two   subroutines   that   are   tailored  to  the  specific  problem,    called 
DIFFUN    and  MATSET,    for  which  examples    are    given   in  the  Appendix.      The 
purpose    of  DIFFUN   is  to    calculate   the  vector  f  with  the   value   of   f 
being  stored  in  DY(l).      The  Jacobian  matrix  is    calculated  by  MATSET 

and  stored  in  the   array  PW.      On   return   from  MATSET,   the   contents   of 

<-)  tw  (  T  \ 
PW(I,J)    are    (    —    T    )    for  I    from  1   to  N   and  J   from  1  to  NY,    and 

a  I  {  J  ) 

,       3DY ( I )    N 
9YL(J  NY) '    f°r        fr°m  1  t0  N    and  J    from  WY+1  to  N* 


Finally,   two  more   subroutines    are   required  that   do  not 
depend  on   the    specific  problem  to  be   solved.      The   first   of  these 
is   MINV,    a  matrix  inversion   subroutine   that  takes    as   input  the 
Jacobian  matrix  in  the   array  PW  and  provides   as   output  the  inverse 
Jacobian  in  the  same   array.      Since   a  number  of  matrix  inversion 
routines  that  are    generally   available  may  be   adapted,   the  particular 
subroutine  used  in  this   investigation  is   not   given   in   the  Appendix. 
The   second  of  these   two   is   MATMUL,    a  short    subroutine   that   multiplies 
the    function  vector   in  DY  by  the   inverse  Jacobian  in  PW  with  the 
resulting  vector  being  stored  in   Fl  on   return.      A  listing  of  this 
subroutine   is   included  in   the  Appendix. 
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3.      NUMERICAL  RESULTS 

In   order  to  have   a   good  basis    for  comparing  the  numerical  results 
from  INTECH  with  the   results   obtained  from  programs  using   other  methods    for 
the  solution   of  nonlinear  equations,   it  was    decided  to   choose  most    of  the 
examples   from  previously  published  and  widely   available   investigations   made 
by  others,    and  to  make    comparisons  with  these  published  results.      This 
avoids   the  possibility  of  imperfectly  implementing  the  various   authors' 
ideas   into   computer  programs,  which   could  lead  to  invalid  conclusions. 
However,   it    does  not  make    it   possible  to  have  the   comparisons   as   straight- 
forward as  would  be   desirable,   in   some    cases,  because   of  different 
approaches    and  different  kinds   and  amounts   of  detail  being  presented  by 
different   authors. 

In   addition  to   making  comparisons  with  the  published  results 
of  others,   it  was    decided  to   also  run  each  problem  on   an   already  available 
library   subroutine  for   finding  the  minimum  of   a  function  based  on  the 
Fletcher-Powell   algorithm   [8]    as    detailed  by  Wells    [9]    and  later  modified 
by  Fletcher    [10].      This   involved  only  making  modifications   to  translate 
the   available  subroutine   from  FORTRAN  IV  to   CALL  OS   FORTRAN  so  that   it 
could  be  used  in  the  timesharing  mode   along  with  INTECH.      This   modified 
Fletcher-Powell  subroutine  is   referred  to   as   FLPMIN  in  the  discussion   of 
results.      The   function  to  be   minimized  by  this   subroutine  was  taken  to 
be  the  sum  of  the   squares    of  the   functions   that   formed  the   system  to  be 
solved  simultaneously  in  each   case. 
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3.1     Description    of  the   Examples 

Example    1 

Function: 


f(x)   = 


h  +      x,    +      x_   -    x.    +   2xnx0   +    3x, 


1  2 


1  +   2x±   -   3x2   +   x±  +      x1x2   -   2x, 


Starting  Vector:      x  =   (-2.057,   -7-503) 
Approximate   Solution:      x=    (3-339,   -2.98*0 

This    example   is    one   of  those   used  by  van   Melle    [l]    in   evaluating 
the    subroutine   DIFMF3  upon  which   INTECH   is   based. 

Example  2 

Function:      same   as  Example   1 

Starting  Vector:      x  =    (0,1) 

Approximate  Solution:      x  =    (-1.533U,   0.06ll2l) 

It  was   noted  that   Example   1  had  another  solution,   and  this 
starting  vector  was  tried  with   the  expectation  of  obtaining  the 
above   solution. 


Example    3 
Function: 


f(x)   = 


2 

Xl   " 


*2 


+  1 


xn    -   cos  (J-  x_) 


1  _x2    *2 

Starting  Vector:      x  =    (1,0) 
Solution:      x  =    (0,1) 

This  example   is   discussed  in   detail  by  Boggs    [U,5].      The 
solution   given   is   the   one  which   can  be   located  from  this  starting  point 
without   crossing  the   lines   of  singularity   for  the  Jacobian. 
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Example   k 

Function:  same  as  Example  3 
Starting  Vector:  x  =  (-1,1) 
Solution:      (-  — ,   1.5) 

This   particular  starting  vector  was  not   discussed  by  Boggs, 
but   it   is  representative   of  the   starting  region  which  he   states  should 
lead  to  the   given  solution    [H,5]. 

Example  5 
Function: 


f(x)   = 


1  *2 

-  [sinU^)    -  —  -  x1] 


t1  "  57>   t 


2x±  ey^ 

e  -  e]  +  —  -  2ex1 


Starting  Vector:      x  =    (O.U,    3) 
Approximate   Solution:      x  =   (0.3,  2.8) 

This   is  the  second  of  the  two  examples    discussed  by  Boggs    [4,5] 
and  it   is    also   given  in   an  earlier   article  by  Broyden    [ll]. 


Example  6 
Fun  ct  i  on : 


f(x)   = 


10xl  2 

Z^T+2x2 


Starting  Vector:   x  =  (3,1) 
Solution:   x  =  (0,0) 

This  is  one  of  the  examples  discussed  by  Powell  [6],   One  of  his 
earlier  methods  did  not  converge  to  the  correct  solution  of  this  problem, 
owing  to  the  singularity  of  the  Jacobian  matrix  at  points  between  the  starting 
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point  and  the  solution.      Powell  indicated  that  his  new  algorithm  was 
developed  in  order  to  provide    convergence   for  examples   such   as  this, 
but   information   on  the  number  of  function  evaluations    required  for 
this  particular  example  is   apparently  not   given   in  his    discussion   of 
numerical  results    [7]. 


Example   7 
Function: 


f(x)   = 


IOjOOOx^x^ 


-   1 


-x^ 


-X1  "2 

e  +   e  ■-   1.0001 


,-5 


Starting  Vector:      x  =   (0,1) 

Approximate   Solution:      x  =    (1.098  x  10"v9   9.106) 

This   example  is   one  indicated  by  Powell    [7]    as  being  difficult, 
primarily  because   of  being  badly  scaled. 


Example   8 
Function: 


f(x)    = 


10( 


-   x. 


■*2        "I' 
1   -\ 

Starting  Vector:      x=    (-1.2,   1.0) 
Solution:      x  =    (l,l) 

This   example  was  used  by  Powell    [7],   and  appeared  earlier  in   an 
article  by  Rosenbrock.    [12].      Reference   is    also  made   to   it  by  Broyden    [ll] , 
Branin    [13],    and  Fletcher  and  Powell    [8]. 


f(x)    = 


Ik 


Example   9 
Function: 

"x1(x1(5  -  xx)    -     2)   +   x2   -  13" 
x1(x1(l  +   x1)    -  lU)   +  x2   -  29 

Starting  Vector:      x  =    (15,   -2) 
Solution:      x  =   (U,5) 

This  problem  is    given   in  the  Appendix  to  Powell's    article    [7] 
as    an  example   of  one  which   does  not    converge  "by  his   method. 

Examples   10,    10a,    11,    11a,   12,    12a,    13,   13a 

Function: 

n 

£      (A.  .    sinx.   +  B.  .    cosx.)    -  E.    =   0,      i=l,2,.  .  .n. 

The  elements   A.  .    and  B.  .    are    calculated  by   a  pseudorandom 
number   generator  intended  to   give   uncorrelated  numbers    from  a  rectangular 

distribution  between   -100   and  +100.      The  numbers  E.    are    calculated  to 

i 

satisfy  the  equations    for  a  particular  vector  x*  =    (x*,    x*,   x*,...x*) 
where   each   of  the   components   was   selected  by  the  pseudorandom  number 
generator   from  a  rectangular   distribution  between   -tt  and  tt  . 
Starting  Vector:      x  =   x*  +   O.ln,  where   the   components  of  n  =    (ru  ,   ru , .  .  .n    ) 
were   selected  by  the  pseudorandom  number   generator  from  a  rectangular 
distribution  between   -tt  and  tt  . 

Solution:      x*  =   ( x* ,    x*,...x*)    as    described  previously   for  this    functi 
For  Examples   10   and  10a,   n=5;    for  Examples   11  and  11a,   n=10; 
for  Examples    12   and  12a,   n=20;    for  Examples   13  and  13a,  n=30.      The  two 
examples    for  each  value   of  n  have   different   coefficients   obtained  by 
starting  the  pseudorandom  number  generator  at   different  points   in  the 


on, 
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two   cases.      Results   of  the   solution   of  examples   of  this   nature  were 
reported  by  Powell    [7]    and  by  Fletcher  and  Powell    [8]. 

3.2      Discussion    of  Results 


The   results    from  computer  runs    of  the  seventeen  examples   using 
the   van   Me  lie   method   (IN  TECH)    and  the   Fletcher-Powell  method   ( FLPMIN) 
are   summarized  in  Table   I    and  Table   II.      Also  shown   for  comparison    are 
the  published  results    from  Boggs    [H,5]   in  Table   I,   and  from  Powell    [7] 
in   Tables   I   and  II    for  those   of  the  seventeen  examples   that  were  used 
by  them.      Table  I   summarizes   the   results    from  the  examples  having  just 
two   equations    in   two   unknowns,   while   Table   II    includes   those   examples 
having  from  five  to  thirty  equations    and  unknowns. 

In  Table   I  the    convergence    criterion   for  INTECH   is   of  the 
same  form  as   that  used  in  van   Me  lie '  s   subroutine   DIFMF3,    and  the  value 
is   chosen   so   that  the    criterion  will  be   more  stringent  than  those 
(each   different)    used  by  Boggs   and  by  Powell  in  their  examples.      For 
FLPMIN  the   criterion  was    chosen  to  be  of  the   same    form  as   that  used  by 
Powell  since   the   sum  of  the   squared  functions   is   needed  for  other 
purposes   in  both  methods,   but  the   value   used  was    designed  to   make  it 
about   as    stringent   as   the   criterion   used  for  INTECH,    assuming  the  error 
is   about   equally  divided  between   the  two  variables.      For  the   examples 
in   Table   II   subroutines   INTECH   and  FLPMIN  were  both   altered  to  use 
exactly  the  same  error   criterion   as  that  used  by  Powell.      Such   changes 
appeared  to  be  necessary  in  order   to  make  valid  comparisons    for  these 
examples,  which   involve   many  more   variables  than  the  ones   in  Table   I. 
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Table   I.      Comparison  of  Results    for 
Examples  with  Two  Unknowns 


EXAMPLE 
NO. 

INTECH 

FLPMIN 

BOGGS 

POWELL 

I 

J 

NFE 

I 

NFE 

I 

NFE 

NFE 

ACC 

1 

16 

1+ 

2k 

35 

105 

2 

11+ 

1+ 

22 

2U 

72 

3 

21 

U 

29 

27 

81 

17 

36 

k 

T 

2 

11 

F 

- 

5 

16 

1+ 

2U 

26 

78 

18 

36 

6 

31 

6 

1+3 

F 

•    - 

7 

ll+6 

26 

198 

F 

- 

22  3 

10-10 

8 

12 

3 

18 

66 

198 

28 

10-6 

9 

102 

13 

12  8 

30 

90 

F 

ID"6 

I   -  Number  of  iterations 
J  -  Number  of  Jacobian  evaluations 
NFE  -  Number  of  function  evaluations    or  equivalent;   a  Jacobian   or 
gradient  evaluation   is   counted  as  n   function   evaluations 
n   -  Number  of  unknowns 


Convergence    Criteria 


INTECH:         Z       |  f X  \    <   10 
i=l 


-6 


i    ?  -1? 

FLPMIN:         Z      (f   )      <    .5   x  10 

i=l 


BOGGS : 


f1!    <   10"5,   i=l,2,...n 


i    ? 
POWELL:         Z      (f   )      <  ACC 

i=l 


IT 


Table   II.      Comparison  of  Results    for 
Examples  with  More   than  Two  Unknowns 


EXAMPLE 

INTECH 

FLPMIN 

POWELL 

NO. 

n 

I 

J 

NFE 

I 

NFE 

NFE 

10 

5 

9 

3 

2k 

33 

198 

11 

10a 

5 

9 

3 

2k 

39 

23U 

12 

11 

10 

13 

3 

1+3 

k6 

506 

19 

11a 

10 

10 

2 

30 

k6 

506 

23 

12 

20 

13 

3 

73 

79 

1659 

36 

12  a 

20 

F 

- 

- 

77 

1617 

36 

13 

30 

F 

- 

- 

107 

3317 

U7 

13a 

30 

19 

3 

109 

113 

3503 

U3 

I   -  Number  of  iterations 
J  -  Number  of  Jacobian  evaluations 
NFE   -  Number  of  function   evaluations  or  equivalent;    a  Jacobian   or 
gradient  evaluation   is    counted  as   n   function  evaluations 
n   -  Number  of  unknowns 


Convergence   Criterion 

For  all  methods   in   the  table : 


n 


(f1)2    <    10"3 


i=l 
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The   subroutine  INTECH   fails   to   converge   for  two   of  the   seventeen 
cases,  Examples   12a  and  13   as   shown  in  Table   II.      It   does    converge, 
however,   for  Examples   12    and  13a,  which   are   of  the  same  nature,   involving 
twenty   and  thirty  variables,    respectively,  but  formed  from  different 
sets   of  random  numbers.      Apparently  the   method  of  determining  the  values 
of  H   and  ALPHA  to  be   used  is  not    completely   satisfactory  in  these  two 
cases  because   the  error,    after  decreasing  for  awhile,    increases   markedly 
for  several  steps   then   gradually  and  continues  to   cycle   in  this   manner. 
It  was   possible  to  obtain   convergence   for  Example  12a  by  making  a  change 
in  the  method  so  that   slightly  smaller  values   of  H  were  used.     This 
change   increased  the  number  of  function   evaluations   required  in   other 
cases,  however,    and  no   further  attempts  were   made  to   optimize  the 
strategy  because   the  primary  purpose   of  the  investigation  was  to 
evaluate   the  method  as   it  had  been   developed  up  to  this  time. 

The  Fletcher-Powell  subroutine,   FLPMIN ,   applied  to  the 
function   formed  by  taking  the  sum  of  the   squares    of  the  original  functions, 
failed  to   converge   to   a  solution  of  the   original  equations    for  Examples 
k,    6,    and  7-      In  the    case  of  Example   k  it   converges  to   a  nonzero  minimum 
of  the   function,  which  is   not   a  solution  of  the  original  equations.      In 
the    cases   of  Examples   6   and  7  it  will  be  necessary   for  the   reader  to 
refer  to  Fletcher's    algorithm   [10  ]    for  the   definitions   of  the  names 
given  in  quotation  marks  in  the   following  discussion.      In  Example   6  the 
reason   for  non- convergence  was  that  the  variable   "gb"    as    calculated  just 
after  the  label   "search   along  s:"   becomes  positive,  which  means   that  the 
matrix  "h"    is  no  longer  positive   definite    as  it   should  be   if  the 
procedure    for  updating  it  were  exact    [8].      In  Example   7  the  problem  is 
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caused  by   continuous    looping  through   the   section   "beginning  with  the   label 
"interpolate:",  without  ever  satisfying  the   criteria  for  leaving  that 
loop.      It   is  possible   that  the   difficulties   encountered  in  these   latter 
two  examples    could  be    circumvented  by   changes   in  the  subroutine,  but  no 
attempt  was    made   to   do   that   in   this    investigation. 

In  the   twelve  examples   in  which  both  INTECH   and  FLPMIN   converge, 
INTECH  requires   fewer  function  evaluations   in   all  except  one   case.      On 
the   average  it   requires   less   than   one-half  of  the    function  evaluations 
required  by   FLPMIN   in   the   examples   of  Table   I,   which   involve   only  two 
variables.      In  the  examples   of  Table   II,  which  involve   more  than  two 
variables  ,   FLPMIN   seems  to   require    on  the   average   about  n-times   as   many 
function  evaluations    as   INTECH  where  n   is   the  number  of  variables   in  the 
example. 

Examples   1   and  2   are  the   same  problem  with  different   starting 
vectors  that  were   expected  to  lead  to   different  solutions.      The  same  is   true 
of  Examples    3   and   k.      For  Examples    1    and  2    INTECH  produced  the   expected  solutions 
in  each   case  while   FLPMIN    converged  to   the  same  solution   (the   one  expected 
for  Example  2)    in  both   cases.      For  Example    3  it    converged  to  the   solution 
expected  for  Example   h,    and  for  Example   U — as   shown   in  Table   I — it    found 
a  nonzero  minimum  which  is  not  a  solution   of  the  original   equations. 
Actually,   it  is  not   surprising  that   such   differences   should  occur  "between 
the   two   methods   since   the    function  being  minimized  in   FLPMIN   is  the  sum 
of  the  squares   of  the   functions   used  by  INTECH,   and  such  a  transformation 
could  be  expected  to  produce   different  regions   of  convergence  to   the 
various    solutions.      In  the  practical   application   of  any  method  to    complex 
nonlinear  equations   that  have  not  been  previously  solved,   the  problem  of 
determining  an  appropriate   starting  vector  for  that   particular  method 
will  exist. 
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Results    from  Example   3   and  Example    5  were   also    given  by  Boggs    [,U,5] 
Table   I   indicates    that    INTECH   requires    from  twenty  to   thirty  percent   fewer 
function   evaluations    than   the  Boggs   methods.      Actually,    the   table   shows 
the  best   results   obtained  by  Boggs  with  each  of  the  two   examples   out   of 
twelve    different   implementations   of  his     basic   method  that  were    applied  to 
each  example.      It   turns    out   that   his  best    algorithm  for  Example   3   did  not 
converge   for  Example    5    and  vice-versa.       Considering   only  those   of  his 
algorithms  which   converged  for  both  of  his   examples,   the  best  one   requires 
more   than   twice   as    many   function   evaluations    as    INTECH   in   each    case. 

Referring  to   the   results   reported  by  Powell    [T]    for  Examples   7, 
8,    and  9?   it   can  be   seen   in  Table   I  that   INTECH   compares    favorably, 
requiring  fewer   function  evaluations    for  Example   7  and  Example    8,    and 
converging  where  Powell's   method  does   not    for  Example  9.      Powell's   method, 
however,   requires    consistently   fewer   function   evaluations   than  INTECH   for 
Examples    10  through  13a  as   shown   in  Table   II.      Comparisons    cannot  be 
completely   straightforward  in  these  examples,    of   course,   because   the 
coefficients   were   generated  by  pseudorandom  number  generators    and  were 
probably  entirely  different   in  the   two    cases. 

After   completing  the    computer   runs   on   all  the   examples,    it  was 
pointed  out   to  the    author  that   an   example    from  Rosenbrock    [12  ] ,   which 
was    discussed  in  a  recent  article  by  Branin    [13],  was   one    for  which 
certain  methods    did  not    converge.      The  reference  was    found  and  the  example 
was    quickly   tried  on  both  INTECH   and  FLPMIN  ,   requiring  twenty-nine    function 
evaluations  with  the   former  and  19 8  with  the   latter.      Subsequently,   the 
author  realized  that  this   example  had  previously  been   done  as  Example    8, 
having  been    found  in  Powell's    article    [7]   in  the   first    instance.      On 
comparing  the  number  of   function  evaluations   required  in  the  two   runs   of 
the   same   problem,   however,    it  was    found  that   Example    8   only   required 


21 


eighteen   function  evaluations  with  INTECH,   although  with  FLPMIN   it 
required  the   same  number,   198,    as    on  the   later  computer  run.      The 
discrepancy  was   resolved  when   it  was  noted  that  this   example   is    one  in 
which  one   of  the  variables,   called  x     in  Example   8,    appears  only 
linearly.      In    doing  Example   8  this    fact   had  "been  properly  noted  and  the 
special  provision   in  INTECH   for  linear  variables  had  been  utilized  in 
implementing  the  example.      When   the  same   example  was   inadvertently 
used  the   second  time,   the   linearity  of  x     was    overlooked  and  both 

variables  were   treated  as   possibly  nonlinear  ones   in  the   implementation 
for  INTECH.      Thus,   the  reduction   in  the  number  of  function  evaluations 
required  in  van  Melle's   method  when   advantage   is   taken  of  the  special 
provision   for  linear  variables  was    documented,   although  this  had  not 
been   planned  as   part    of  the   investigation. 

Example  9,   the   one   for  which  Powell's   method  did  not    converge 
although  success  was    obtained  with  both  INTECH   and  FLPMIN,  was    given   as 
an  example   in   an  earlier  article  by  Freudenstein   and  Roth    [l^]. 
Convergence  was   obtained  by  their  method,   although  the  total  number  of 
function  evaluations   required  was  not    given.      Using  a  parameter 
perturbation  technique,  ten   steps  were   required,   each  step   involving  the 
solution  by  Newton's   method  of  a  nonlinear  system  having  different  parameters 
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k.      CONCLUSIONS 

The  primary   conclusion   is   that  the   integration  technique 
developed  by  van   Melle    for  the   solution   of  nonlinear  equations   has 
considerable  merit.      As   implemented  here    in  the  subroutine  INTECH  ,   it 
was  tried  on   a  number  of  different  examples.      Most   of  these  had  been 
identified   in  previous  publications    as    examples  which  were   known   to 
be   troublesome   to   solve,    at   least  by   some   methods.      Every   example   that 
was   tried  on   INTECH   is    reported  in   Section    3.      Under  these    conditions, 
the  fact   that   convergence  was    obtained  in   fifteen   of  the   seventeen 
cases   tried,    and  that   the  number  of   function   evaluations   was   less   than 
for   the   other  methods    compared  in    about   half  the    cases  ,    gives    support 
for  this    conclusion. 

Another   conclusion,   based  on  thinking  about   the  situation 
rather  than   on  experimental   data,    is  that   it   is  probably   feasible  to 
further  improve   the  implementation   of  this   method.      The  procedure  used 
to  establish   the   step   size   H   and  the  parameter  ALPHA   for  the  next 
iteration  was    determined  empirically;    and  since    it  works  well  in  many 
cases  ,   it   should  be   retained  in  basically  the  same    form  unless 
extensive    additional  testing  proves   some   different   approach  to  be  better. 
It   is   suggested  that  by   studying  in   detail   all  examples    for  which   the 
method  is   found  not   to    converge,   or  to    converge  slowly,   it  may  be 
possible   to  supply   additional  logic  to   recognize   in  the   early  stages 
similar   cases    and  take   more   intelligent   action   for  each  of  such 
specific  situations   studied.      Examples   in  which  these  particularly 
troublesome  situations    did  not    develop  would  be  handled  by  the 
present  procedure. 
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The  empirical  procedure   for  deciding  when  to   re-evaluate  the 
Jacobian   and  its   inverse   seems  to  he  reasonably  effective   in   reducing 
the  number  of  function  evaluations   in  INTECH ,   as   indicated  particularly 
by   comparing  with  the   results    from  FLPMIN    in  Table   II,  which   includes 
all  the  examples   involving  five   or  more  variables.      However,   it    is 
possible   that   matrix  updating  procedures   such   as   those   described  by 
Broyden    [15]    or  the   one  used  by  Powell    [6,7]  would  prove  to  be  even 
more   economical  of   computer  time,  especially  for  systems  with  many 
unknowns.      This  possibility  is   supported  by  the   fact   that   in  the 
examples   of  Table   II,   INTECH   consistently  requires   more   function 
evaluations  than  Powell's   method,    and  the  number  of  these  used  for 
the  Jacobian  was    greater  than  the  number  used  for  iterations.      It   is, 
therefore,    concluded  that  the   use   of  matrix  updating  procedures   in 
combination  with  the  basic  approach  of  van  Melle   may  be  well 
worth   investigating. 
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) , YLSV<2) ,OEL 


M*  EXAM°LE     g 

PFAL*S    Y(2),nY(2),Fl( 2),YL(2) ,SAVE(2, 
PEAL*3    YD<2) ,PWQRK(2) 

DIMENSION    PW (?,?) 
NCX  =  8 

Y<1 ) =-1.2 DO 
YL( I )=1.00 
M=2 
NL  =  1 

DEL=1.0-6 
NIY  =  N-NL 
NSFMD=110 
WRITE (6,1  J     MEX 
1     FOPMAT     (i-PX&vp|_c    Nn.«,ic///) 

CALL     INJTFCH{Y,YL,nY,PW,DPL,Pl,YO,  S  A  VE  ,  YL  S  V  ,  P  WOR  K  ,  % 

NtNYfNL tNSEND) 
END 
SURROt)TI  NF    HI  F«=UN  (  D  Y,  Y  ,  YL  ,  N,  NY  t  NL  ) 
PEAL*3     DY(M), Y(NY) ,YL(NL) 

DY(1  )=10.D0*{ YL(1 )-Y<l  )**2) 
DY( 2)=1.D0-Y( 1) 
Q  PTLPN 
END 
SUBROUTINE     MAT  SCT{  PW,Y,  YL  t  N,NY,NL) 
REAI  *=>     Y(MY)  ,  YL  (  NL  ) 
DI  MENS  ION    PW(N,N) 

PW( 1  ,1 ) =-20.D0*Y(l ) 
pwi ! ,2) = !0.no 

PW( 2, 1  )=-l .00 
PW( 2,2) =0.00 

R  E  T  UP  N 

END 
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SUBPDUTINF     INTECHi  YtYUDY,  PWtDEL  tpl  tYD,SAVE,  YLSV% 
, PWO°K,N»NY, NL.NSEND) 

REAL*B    DY(N),F1(N),Y(NY),YL(NL),YD(NY),SAVE(2,NY),YLSV(NL)?; 
, PrfORM N) ,0MAX1,DMIN1 ,  DEL, ALPHA, HfR,DET,SStS ,HOLD ,DD , D  ,  HH 

DIMENSION    PW(N,N) 

11  ************************************************************ 
•«    PARAMETFR    LIST 

"  Y        V=rTRR     OF     M1NLINEAP     VARIABLES     (INPUT). 

»  YL        VECTOR     OF    LINEAR    VARIABLES     (INPIJT). 


ii  ny       VFCT0»    Oc    FUNCTION    VALUES,    CALCULATED 

"  RY     TH^     SUBROUTINE    DIF^UN. 

"  PW       JACOB  I  AN    MATRIX    CALCULATED    BY    SUBROUTINE    MATSET 

"  ns     thc     INVERSE     JACORIAN     AFTER     CALLING    M[NV. 

»  DFL        FRROP     C.p  ITE"  IOM.SUOPLIEO    BY    CALLING     PROGRAM. 

•»  F!        VECTOR     EDUAL     TO    PRODUCT     OF     PW     AND     DY . 

ii  yd       THC    VECT-^P    D^    CHANGES    TQ    BE    APPLIED    TO    VECTOR    Y. 

"        SAVF       TWD-POW    VECTDD    USED    TO    ST^PF     VALUES    OF    Y    AND    YD 
m  rc^ORE    niexT    ITc5^j^fj,     i nj    C4sc    RF-START    Nr  =  DED. 

"       YLSV       V  =  CTnc?     USFD    TO    STOPE     VALUES    Oc    YL    BEcOPP 
"  NEXT     ITERATION,     IN    CASE    PE-STAPT    NEEDED. 

••  PwnPk  USED    "INLY    cOR     STORAGE    SpACE    WITHIN    SUBROUTINE    minv. 

*  N  TUP    TDTAL    NUMBER     OF    VARIABLES     (INPUT). 

"  NY  T"c    NUMBER  OF  N1NLINPAP  VA°IABLES  (INPUT). 

"  NL  TUP  NUMBCR  Oc  LINEAR  VARIABLES  (INPUT). 

'•  NSpND  THE  maximum  NUMBER  DF  ITERATIONS  DESIRED  (INPUT). 

!•*****************************#-****************************** 
l»  ********** 

"     INITIALISE    VARIABLES 
11 

\JDCAT!   =1 
JACOR=0 

p  =1 

N  7  -  M  *  N 
NJ=NY*2 
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I F(MJ.LT.!0)    NJ  =  10 

NJJ=NJ*2/3 

ALPHA=0 

H=l 

MS=0 

NW=0 
i»  *******  *** 

'»    PQINT    HEADINGS 

M 

MRITE(6,1)    N,NL,OEL 

1  FORMAT(»iN    =»,I4,«  NL    =»,I^,«  DEL     ='»D10.2) 
WRI  TE(6,2) 

2  FOR^ATCO    MS    NW       ALPHA  H  ',  8X,  •  ERROR  '  t  6X,  % 

»  Y(  J)     AND    YL(* ) •// ) 

II  **  **  ****** 

"     INITIAL    EVALUATION    nF    FUNCTION    (DY). 
if 

NS  =  NS+! 

CALL     DIFRJN(  DY,Y,YL,N,  NY,  ML  ) 
ii  ********** 

"     INITIAL     EVALUATION    qf     JAC06IAN     (PW). 
ii 

CALL    MATSFT(PW,YtYL, Nf NY, NL I 

************ 

»•      INITIAL     cVAH)ATinN|    If     INV^S^     JACO^IAN     (  PW  )  . 
ii 

NW=NW+l 

CALL    MTNV({JW»!DET,N,NZ,Fl,PWORK) 

«  *  *  ******** 

"     PRODUCT    °W*ny     IN     V^CTPR     Fl. 
•• 

C  A  L 1     M\TMUL(PW,DY»F1, N J 
•i  *******  *  ** 

"    NORM    STORED    IN    SS    FOP     STATING    VALUES. 


03    ?    J  =  !  t  N 
3    SS=SS+DARS ( ^Y ( J  >  ) 
ii********** 

••     pcinjt     RFSULTS    OBTAINED     4ITH     STARTING     VECTDR. 
ii 

WP ITC( s,i ? )    NS,NW, ALPHA,H,SS,(Y( J) , J=1,NY) 
Ic(Nt.  .GT.fl)     WRITE!  6  t1-^)     (  YL(  J)  ,  J  =  1,ML) 

II********** 

ii  PREPARE  FOR  THF  INITIAL  PREDICTOR  STEP  FOR  Y;  DO  THE 

"  INITIAL  CORRECTOR  STFP  fhr  YL  . 

•i 

DO  4  J=! , NY 

U     Y^(  J  )=-M*ci  (  j  ) 

T  c  <  N  L  .  L  r  .  0  )  GO  TO  6 
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oo    5    J=1,NL 

5  YL( J)=YL( J)-F1 ( J+NY) 

6  CONTINUE 

********** 

SAVE     PRESENT    VECTORS    Y    AND    YL. 

DP    7    J=1,NY 
SAVEd  ,  J»  =  Y(  J  ) 

7  SAVC(2, J)=YD( J) 
IC(NL  .Lr.O)    G^    TO    9 
o  0    fl     J  =  !  f  N  L 

R     YLSVf J) =YL( J) 

9    CONTINUE 

********** 

INOIC4TF  NEWTON'S  METHOD  IS  IN  USE  BY  NEWTON  =  1. 

NEWTQN=! 
********** 

BEGIN    MAIN     ITERATION    LOOP    BY    COMPLETING    THE    PREDICTION    OF    Y. 

10  DO    I  1     J  =  l ,NY 

11  Y( J)=Y( J)+YD< J) 

**  ******** 

FVALUATE  THF  FUNCTION  (DY). 

N  S  =  N  S  + 1 

CALL  DIP-UN(DY,YfYL,NTNY ,NL ) 
********** 

CALCULATE     THE     VALUE     OF    T  HP     Nn<5M     (S). 


S=0 
DO    1  -> 


uu  J  =  1  ,N 

12    S=S*0ABS<DY(  JH 

********** 


********* 

PRINT  RESULTS  FOLLOWING  PREDICTED  STEP,  THEN  RETURN  IF 
NORM  IS  LCSS  THAN  QR  cOUAL  T0  DCL. 


WRI  T>-(^,  !3  )  NS,NW, ALPHA, H,  S,(Y(Jl,J=l,NY) 
13  FORMAT  (14, I3vF6«2«2fUl«2t2Dl5.5/(5Xv4D15«5)) 

IF(NL.ST.O)  WRITE! 6, 14)  (  YL(  J)  ,  J--!,NL) 
] U    F 00  MAT  (^X,^D!5.5) 

IF(S.LF.DFL)  RETURN 

**  ******** 

CALCULATE  PREP  (PATIO  OF  NOpM/ °R c V I OUS  NORM). 

PR  ED=S/SS 
********** 

BEGIN  PROCESS  OF  DETERMINING  VALUES  op  ALPHA  AND  P 
BY  WHICH  H  WILL  BE  MUL.TIPLIEO)  e?R  N-XT  I  Tc RATION. 


( FACTOR 
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GO 

TO 

35 

on 

Tl 

IS 

GO 

TO 

30 

IF ( NEWTON. EO. 1) 
IF< PREO.LT.100) 
IF(NOFAIL.Eg.l  ) 

!5  CONTINUE 

IF (PPE0.LT..9R)  GO  TO  16 

ALPHA=l 

PaDMiMl ( 1. 300, 0. 600/H ) 

GO  to  17 

16  D=1.7DT-.e5D0*H+,! 5D0/H 
ALPHA=ALPHA*0.8 

"  SAVF  PRESENT  VECTORS  AND  VALUE 
it 

17  DO  19  J=1,NY 
<^AVE( 1, J)=Y( J) 

IB  SAVE(?, J)=YD( J) 

IF(NL.LF.O)  GO  T0  2  0 

0^  19  J=1,NL 
19  YLSV( J>=YL( J) 
2  0  CONTINUE 

HOLD=H 

•«  ********** 

ri  DcTcpMlNr  WHETHEP  J  AC  OR  IAN  SHOULD 
IF( SS.GF.l.OOJ  GO  TQ  21 
!»=(  JACTR.LT.NJJ  )  00  T3  ?-> 
21  CONTINUE 

JAC  3B  =  JAC0R-1 

IF(  JAC">B.GT.0)  GO  TH  ?3 

************ 

"  EVALUATE  THC  JACOBIAN  (OW), 


OF  H  BEFORE  NEXT  ITERATION. 


BE  RE-EVALUATEO. 


22  CALL  MATS-T(PW,Y,YL,N,NY,NL) 

M w=N W+l 

»♦  *  *  *  *  *  *  *  *  A  * 

"  EVALUATE  the  INVERSE  QF  JACORIAN 


(PW)  . 


CALL     MINV(PWfOET,NtN/,Fl,PwnPK) 
J  A  C  0  B  =  N  J 


f  0  D  m   v  e  c  t  o  p 


THE     PonrjuCT     PW*OYS     THEN    DO    THE 


CORRECT^     STrP     ^OP    Y. 


23     CALL     M*TM|jL(  ^,nY,cl,  N) 
IF(  ALPHA. LT. 0.01  )     GO    Tn 
DD=1  .D0/(  !  *H*ALPHA  > 
qo     -»4    J  =  l  ,  NY 
n=( cl (   I)*H*YD(J) )*D0 
Y(    I  )=Y( J}- ALPHA*D 

?4    YD{ J)=P*(YO(J )-Q) 

GO   t n    ->  7 
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25  HH=-H*R 

nn    26    J=l ,NY 

26  Y0<  J)=HH*FK  J) 

27  CONTINUc 
IP(NL.LE.O)     GO 

"    DO    THE    CORRECTOR 


TO    2o 

STEP    FDR    THF    LINEAR     VARIABLES     IN    YL. 


on    2  8    J  =  1  ,NL 
2B    YL(  J)  =  YL( J  )-Fl  (J+NY) 

29    SS=S 

NOFAIL-1 

Ic( NS.GE.NSENO) 

"    END    Dc    mmm 


GO    TO    3  7 
nation    LOOP. 


GO    TO    10 
30    CONTINUE 

*•    CHANGE     ALpHAf     R, 

ALPHA=! .0 
IP(JACOB.LT.NJJ 
Q  =  r>M  AVI    (  0.500,  0 

H=HOLD*P 

"  PREVIOUS  VFfTPRS 
"  RE-STARTING  LAST 

OP  3  2  J=l ,NY 
Y< J)=SAVE(1 ,  J) 
Y0( J )=9AVF( 2, J) 
IF(NL-Lc.O)  GO 
on  ^3  |=t,NL 
YL( J)=YLSV< J ) 
FONT  I  NUE 
MOFAIL=0 
GO  TO  !  0 
CONTINUE 

If********** 

"  IF  (  N F W TH N s ]  A N n 

"    C0PRCCTI0N    PROCES 

IF(P»E0,LT.0.95 

ii  *#  *^ * **  a* * 

"  BUT  IP  PRED>=.9  5 
"  NEWTON,  J AC 01;  TH 
'•  ALTFP  VALUE  n  F  NO 


*1 

3-> 


^3 
3^ 


•*«; 


ANO    H    C0R     THE    CASE     (PRFDMOO    AND    NOFAIL=l) 


)  JACOB=0 

.7D0/HOLn  ) 


RESTORED    AND    NOFAIL     SrT     TO     ZERO    BPF3RF 
STcp,     BECAUSE    nF    LARG-     INCREASE     IN    ERROR 


T0    3^ 


PRFD<.:>E)  GO  BACK  TO  00  THF  NORMAL 
S,  REMAINING  IN  THE  NEWTON  METHOD. 
)  GO  TO  17 

PRIMT     A     MFSSAG-,     CHANGE     ALPHA,     H,     R, 
EN    GO    TO    tFSTORE    PREVIOUS    VECTORS    AND 
FAIL    BEFORE    RE-STARTjsjG    LAST     ITERATION, 


WRI  T=(6, ->6)S,  (  Y<  Jl  ,  J  =  l  ,NY) 
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3  6  poo 

IH 

H=. 

ALP 
R=. 
JAC 
MP  W 
GO 
37  COM 
n  ****** 

"  IF  TH 

"  PRTMT 

W&I 


3  8 
39 
40 


FOR 
WRI 
POR 
WRI 
FOR 
RET 
EMO 


MAT   ( 

NL.OT 

01 

HA  =  1 
01 
0R  =  0 

TOM  =  0 

m  3i 

TIMUF 
**** 

F  MUM 

THP 
TE(6, 
MAT(  • 
TE<5, 

M  \  T  (  • 

TE(6, 
MAT  ( 

UPM 


•  NEWTON  FAILED: • ,9X,011.2, 201 5 . 5 / < 5X ,4015 . 5 ) ) 
.0)  WRITE(6,14)  (YL(J) ,J=1 ,NL) 


BEc?  nc  FUMCTIOM  EVALUATIONS  HAS 

VECTORS  DY,  Yf  YL  ,  !\N0  F!'t  THEN 

3«)     (0Y{ J) ,J=1 ,N) 

-DY: •/( •**•  ,4014.6) ) 

^Q)     I Y( J), J=T ,NY) 

-Y:  •/(•**  '  ,4016. 9.  )  ) 

40)     ( J,F1( J),J=1,N) 

»-F! :«/31 I5,«**« ,014.6) ) 


EXCEEDEO    NSFND, 
3FTIJ9N. 


M  *****  * 
H  **  **** 
•f 

SUBP0UT 

REAL+3 
01  MEMS 

no 

SljM 

on 

sum 

FT  ( 
RFT 

en 


40 
50 


************************* ********* ********************* 
*  *  ***************************************************** 

rfsJE     MATMUL(PW,DY  ,cl  ,N) 

0Y( M) ,F1( M) , SUM 
ION!     PW(N,N) 
*0    J-l ,M 
=  0 

40    K=l , M 

=  SIIM+PW(  J  ,K)  *0Y(  K) 
J  )  =  S!JM 

URN 
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